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We introduce an extension to the ladder truncated Bethe-Salpeter equation for mesons and the 
rainbow truncated quark Schwinger-Dyson equations which includes quark-loop corrections to the 
gluon propagator. This truncation scheme obeys the axialvector Ward-Takahashi identity relating 
the quark self-energy and the Bethe-Salpeter kernel. Two different approximations to the Yang-Mills 
sector are used as input: the first is a sophisticated truncation of the full Yang-Mills Schwinger- 
Dyson equations, the second is a phenomenologically motivated form. We find that the spectra and 
Oh' decay constants of pseudoscalar and vector mesons are overall described well for either approach. 

Meson mass results for charge eigenstate vector and pseudoscalar meson masses are compared to 
lattice data. The effects of unquenching the system are small but not negligible. 

<N 

j>. ■ I. INTRODUCTION 

m 

One of the longstanding problems of the theory of the strong interaction, Quantum Chromodynamics [QCD], is the 
description of hadronic bound states and resonances from first principles. Ironically it is perhaps the most interesting 
aspect of QCD, the phenomenon of confinement, that greatly hinders the development of our understanding of the 
experimental hadron spectrum: the individual components of these hadrons, dressed quarks and gluons, do not appear 
as physically observable states in nature. Thus, in order to study the internal structure of hadrons one has to resort 
to theoretical frameworks that connect the fundamental fields of the theory with phenomenology while preserving the 
basic low energy properties of QCD: confinement and dynamical chiral symmetry breaking. 

Lattice simulations are not entirely satisfactory in this respect. They provide values for the global properties of 
Oh. hadrons (masses, decay widths etc.,) which, when compared with experiment, reassure us that QCD is indeed the 
theory of strong interactions. But they may not be capable to provide enough information on the internal structure 
of these hadrons, which is vital for our understanding of the dynamical aspects of low energy QCD. An alternative 
field theoretical and relativistic method which is well suited to provide this information is the Schwinger-Dyson and 
Bethe-Salpeter formalism (see, for example, 0, 0, for a contemporary introduction to the topic) . In principle 
this approach allows one to derive meson properties directly from the fundamental building blocks of the field theory, 
the Green's functions. 

From a technical point of view, lattice simulations and the Green's functions approach are complementary to each 
other in several respects. Lattice calculations are ab initio and thus contain all effects from quantum fluctuations. 
However, they are limited to a comparably small momentum range and suffer from finite volume effects in the infrared. 
Furthermore the implementation of realistic quark masses turns out to be extremely involved from the viewpoint of 
computational times and is probably not feasible in the near future. Schwinger-Dyson equations [SDEs] can be 
solved analytically in the infrared and provide numerical solutions for a large momentum range. Together with the 
corresponding Bethe-Salpeter framework all aspects of chiral symmetry are respected such that the properties and 
effects of light quarks can be determined systematically with reasonable effort. On the other hand, in order to obtain 
closed systems of equations one has to introduce suitable approximations of higher order n-point functions. These 
truncations have to be controlled by comparing results from the SDE-approach to corresponding lattice data. In turn, 
the Schwinger-Dyson approach provides important information in momentum and quark mass regions that are not 
(yet) accessible on the lattice. 
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One of the most important recent developments in the Bethe-Salpeter formalism is the realisation that the chiral 
symmetry breaking pattern of the pseudoscalar meson sector is governed by the relationship between the quark prop- 
agator SDE and the Bethe-Salpeter kernel supplied by the axialvector Ward-Takahashi identity Q ■ This observation 
has led to the development of phcnomenologically motivated models to describe the light pseudoscalar and vector 
meson properties [3 |7j. These models are constructed using the simplest truncation scheme (the rainbow/ladder 
truncation) known to satisfy the axialvector Ward-Takahashi identity. The Yang-Mills sector of QCD is replaced by 
an ansatz for the effective interaction between the quarks which is sufficiently strong to induce spontaneous chiral 
symmetry breaking and associated dynamical quark mass generation. This provided a robust starting point to deter- 
mine various aspects of the light pseudoscalar and vector mesons like charge form factors HIH, electroweak and strong 
decays 0, and electroweak transitions 0, IT3. Ex plor atory studies go beyond the rainbow-ladder truncation by 
including corrections to the quark-gluon vertex [14. Il5lll6j . 

Despite their phenomcnological success, these models are unsatisfactory in an important respect. It is not clear, 
to what extent the structure of the Yang-Mills sector of QCD is represented by the ansatze used for the effective 
quark-gluon interaction. Certainly, if one strives for a detailed understanding of the internal structure of mesons it 
is vital to investigate (rather than to model) the effects of the Yang-Mills sector on experimental observables much 
more explicitly. This is one of the two central aims of the present work. Based on previous studies 0, Il7| of the pure 
Yang-Mills sector, results for the quenched and unquenched ghost, gluon and quark propagators from their coupled 
set of Schwinger-Dyson equations (SDEs) have been obtained in Ref. 01 ■ Here we extend this setup to include the 
corresponding Bethe-Salpeter equations and determine light pseudoscalar and vector meson observables. The effective 
interaction of the quark and gluons is hereby explicitly resolved in terms of the dressed gluon propagator and the 
dressed quark-gluon vertex. 

The second aim of this work is to quantify a class of unquenching effects in meson observables. In general, 'un- 
quenching' means the inclusion of quark loops to the various correlation functions that enter the expressions for the 
quark self-energy and Bethe-Salpeter kernel. A previous study of unquenching focused on the mechanism for meson 
decay 0| where a single quark loop could be incorporated into the truncation scheme, retaining the crucial property 
of obeying the axialvector Ward-Takahashi identity. However, there exists another viable improvement to the basic 
Bethe-Salpeter truncation scheme and this involves the unquenching of the input Yang-Mills interaction form with 
an infinite series of increasing numbers of quark loops. This is the second main issue of the present work. 

The paper is organised as follows: In section^we detail our treatment of the ghost, gluon and quark propagators. 
We investigate unquenching effects due to quark loops in these quantities in two different setups. On the one hand we 
explicitly solve the coupled set of (unquenched) Schwinger-Dyson equations for the ghost, gluon and quark propagators 
employing ansatze for the various vertex functions. On the other hand we investigate the same effects employing a 
model for the quenched Yang-Mills sector of QCD, which we effectively unquench by adding a quark-loop to the model 
gluon propagator. The results in both setups are used as input for our calculation of pseudoscalar and vector meson 
properties, detailed in section ITTT1 Numerical results and conclusions are reported in sections IIVI and Ivl Various 
technical aspects are discussed in the appendices. Throughout the various stages of our calculation we will compare 
our results from the SDE/BSE approach with those from recent lattice simulations. Finally, it is worth noting that 
whenever we refer to unquenching in the present work we only refer to unquenching effects due to the presence of 
quark loops in the gluon-SDE. Certainly, as already mentioned above, unquenching the full theory involves other 
important effects like meson decays which are beyond the scope of this study. 

II. UNQUENCHED PROPAGATORS 

To start with, we describe in some detail the two calculational schemes that we use to determine the effects of 
quark loops on the propagators of QCD. The first scheme employs the full set of coupled unquenched SDEs for 
the ghost, gluon and quark propagators together with appropriate ansatze for the ghost-gluon, the three-gluon and 
the quark-gluon vertex. We use a slightly modified truncation compared to the one used in Ref. [Tg- The second 
scheme employs a model for the quark-gluon-interaction, which has been shown in Ref. to be well suited for 
phenomcnological investigations. In both schemes we take unquenching effects in the gluon polarisation explicitly 
into account. 



A. Propagators from the quark, gluon and ghost Schwinger-Dyson equations 

For the purpose of our investigation we work in Landau gauge QCD, with the gauge fixing performed by the familiar 
Faddeev-Popov method. The resulting determinant is conveniently represented by a pair of Grassmann-fields, the 
Faddeev-Popov ghosts. The two-point functions of Landau gauge QCD, the ghost, gluon and quark propagators can 
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be generically written as 1 



ur/ip) = - ( ^s a \ (i) 



p 



D*(p) = (^-^^V, (2) 

1 _ Z){p 2 ) 

Slijp) = =ip^)+^P) = -^ + Mi(p2) ■ (3) 

Here G(p 2 ) denotes the ghost dressing function, Zip 2 ) the gluon dressing function, and Ai(p 2 ) and Bi{p 2 ) are 
the vector and scalar quark self energy of the quark flavors i G {w, d, s}. In a somewhat different notation, 
Mi{p 2 ) = Bi(p 2 ) j Ai(p 2 ) is the quark mass function, whereas Z^(p 2 ) = 1/Ai(p 2 ) is called the quark wave func- 
tion. The corresponding bare tree level propagators, D G (p), D° (p) and Sf(p) are given by Z = G = Zf = 1 and 
Mi(p 2 ) = mi, where the mi are the respective renormalised quark masses in the Lagrangian of QCD. The color factors 
of objects in the adjoint representation are given explicitly, whereas the color factors of objects in the fundamental 
representation are treated implicitly. 

The propagators QJ, and |J5J satisfy the coupled set of renormalised Schwinger-Dyson -equations, 

[D af (p%l = Z 3 [D Q ^{p)]-l-Z x J dkf^ bc (p,k)D b G d (p-k)f d ^(k,p)D c J(k) 

-Zx \Jdk t°^ c (p, k) D b p d p ,( P - k) r d ;l a ,(k, P )DZ,(k) 

- Z ^\ E / dk Tri^SiWrXhrfSifr-k)} + (...), (4) 

■i— u,d,s 

[D a G f (p)]-i = Z^D^ipT^-Z, J dkf^ hc (p,k)Dli(p^k)ft e f(k,p)D^(k), (5) 

S- 1 (p) = Z 2 [Sf (p)]- 1 -Z 1F J dk (p, k) S, t (k) Tl (fc, p) Df v (p-k), (6) 

which are given diagrammatically in Fig. ^ Here dk is an abbreviation of d fc/(27r) 4 . The sum in the quark-loop 
runs over the light quark flavors and the trace is over color and Dirac indices. The ghost-gluon vertex is denoted by 
ff' c (k,p), the three-gluon vertex is given by Tp bc a (k,p) and T®(k,p) is the quark-gluon vertex. The corresponding bare 
quantities have an additional superscript 0. The ellipsis in the gluon-SDE represents two-loop diagrams containing the 
four-gluon- vertex. The gluon tadpole diagram vanishes in the process of renormalisation and is therefore omitted from 
the start. The SDEs 10} , JSJ) and JHJ are renormalised and therefore finite and independent of the ultraviolet regulator 
A. The definitions of the renormalisation constants Zi are given in appendix ^ together with the corresponding 
relations between unrenormalised and renormalised dressing functions. 

The quenched system of ghost- and gluon-SDEs, Eqs. JSJ and SJ without the quar k loop, and the quark-SDE © 
have been investigated separately in a series of publications |2l I lL l,5l 1 1 71 |2(I l2ll l22l . 1 2.1 12 ll l2ol |2( il 1271 l'28l |29| . Results 
for the unquenched system, including the back-coupling of the quarks onto the Yang-Mills sector, have been reported 
in [Tq|. In the following we shortly summarise some results, which are important for the present work. 

The infrared behaviour of the ghost and gluon propagators can be determined self-consistently from the ghost-SDE 
and the ghost-loop in the gluon-SDE alone. This so called ghost-dominance of Yang-Mills theory in the infrared has 
been conjectured in Ref. [2(|, lead to the formulation of an infrared effective theory based solely on the Faddeev- 
Popov determina nt E2al and has recently been shown to hold in a self-consistent scenario of solutions for general 
n-point functions [26j. In Landau gauge, such an analysis is based on the factorisation property of the ghost-gluon 
vertex, which entails that in the presence of one external scale this vertex cannot develop any nontrivial dressing (i.e. 
singularities or zeros) in the infrared (22,|3ll3l- Based on this observation one finds simple power laws, 

Z(p 2 ) ~ (p 2 ) 2K , G( P 2 ) ~ ( P 2 )- K , (7) 

for the gluon and ghost dressing function with exponents related to each other. The relations J7J) can be determined 
from the ghost-SDE alone and are independent of any approximation of the ghost-gluon vertex, which is in agreement 



We work in Euclidean space with scalar product a-b = S^a^b^. The Dirac matrices are Hermitian and obey {7^,7^} = 2<5 M „. 
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FIG. 1: A diagrammatical representation of the coupled system of ghost, gluon and quark Schwinger-Dyson equations. Filled 
blobs denote dressed propagators and empty circles denote dressed proper vertex functions. 



with the factorisation property. The exponent k is an irrational number and depends only slightly on the dressing 
of the ghost-gluon vertex |24|. With a bare vertex one obtains k = (93 - %/1201)/98 w 0.595 [231131 This result 
has been confirmed inde pen dently in studies of the exact renormalisation group equation |3lj |. Lattice simulations 
indicate that k « 0.5 |32l l33t 133. l35l l3r| . see however [37j for a discussion on systematic effects on k due to a compact 
manifold. 

Probably the most interesting consequence of the power laws @ is the corresponding fixed point of the running 
coupling a(p 2 ) in the infrared. A nonperturbative definition of this coupling can be derived from the ghost-gluon- 
vertex 20] , 

"(P 2 ) = |W)Z( P 2 ). (8) 

No vertex function appears in this definition; a fact that can be traced back to the ultraviolet finiteness of the ghost- 
gluon vertex in Landau gauge. Note that the right hand side is a renormalisation group invariant, i.e. a(p 2 ) does not 
depend on the renormalisation point. From the power laws (JJJ one finds that the coupling has a fixed point in the 
infrared, given by |24| 

= iL r(3-2 K )r(3 + «)r(i + «) w 

V ; 6N C T 2 (2-k)T(2k) y ' 

for the gauge group SU(3) and a bare ghost-gluon vertex. Interestingly enough, both, the values of the exponent k 
and the fixed point a(0) can be shown to be independent of the gauge parameter in a class of gauges that interpolate 
between Landau and Coulomb gauge |3^ |. 

Numerical solutions for the ghost and gluon propagator for general momenta have been obtained in refs. [Til Hsj. 
In the infrared these solutions agree with the infrared asymptotic behaviour of Eq.J7|), whereas in the ultraviolet 
they reproduce corresponding results from resummed perturbation theory. The employed approximation scheme for 
the ghost-gluon- vertex, the three-gluon- vertex and the four-gluon- vertex is detailed in ^tJ- F° r the convenience of 
the reader we summarise the basic assumptions of this truncation in appendix [5] Here we discuss our truncation 
of the quark-gluon-vertex in more detail, as this vertex provides the crucial link between the Yang-Mills SDEs, the 
quark-SDE and the Bethe-Salpeter equation describing mesons as bound states of quarks and antiquarks. 
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The structure of the quark-gluon-vertex has been the focus of some interest in the last years. It is well known, that 
the leading 7 M -part of the vertex provides sufficient structure to the quark-gluon interaction to describe a range of 
properties of the pseudoscalar and vector mesons well (see |4( and references therein) . On the other hand it is also 
well known by now, that further tensor structure in the vertex, in particular a scalar part, is mandatory to describe 
other meson channels. Such corrections to the vertex have been investigated on a semiperturbative level in several 
works fiH ITHl flfi! Isfj, l4?i Elj and comparisons with first lattice results [32, are encouraging (although the error 
bars from such simulations are still quite large). Unfortunately, the numerical complexity of the present investigation 
is such that we are not able to include such corrections. Instead, we will adopt the strategy of Ref. 18] and employ an 
ansatz for the nonperurbative dressing of the 7 M -part of the vertex such that important constraints like multiplicative 
renormalisability and a correct perturbative limit of the vertex are satisfied. 

With k denoting the gluon momentum and p, q the two respective quark momenta, two possible ansatze for the 
vector part of the quark-gluon-vertex are 

a 2 (k 2 \ 

rt(q,p) = t gT a lv ^->A(k 2 ), 

KM=i 9 T^^^^l±^. (10) 

The presence of the vector quark self energy A(k 2 ) in these expressions is enforced by the Slavnov- Taylor identity for 
the quark-gluon vertex (given in appendix^), which also accounts for one of the two ghost factors. The second ghost 
factor leads to a perturbative behaviour of the vertex, which produces the correct perturbative anomalous dimensions 
of the ghost, gluon and quark propagators in the ultraviolet momentum regime, provided the momentum arguments 
of the dressing functions are chosen in an appropriate way. As discussed in detail in Ref. 0], the first version of 
arguments, is necessary in the quark-SDE, whereas the second version has to be chosen in the quark loop of the gluon- 
SDE to reproduce the respective correct perturbative limits of these equations. Via the renormalisation properties of 
the contributing dressing functions and the vertex renormalisation constant Zip this ansatz has the correct dependence 
on the renormalisation scale fi 2 in agreement with the renormalisation group. Furthermore note that both versions 
of the vertex have the correct properties under charge conjugation. The price to be paid with this construction is 
its cutoff-dependence via Zip. Thus these ansatze do not have all properties of well defined renormalised Green's 
functions, but have to be regarded as effective constructions, which only make sense in connection with the quark- 
and gluon-SDE. If inserted into these equations, the ansatze lead to properly renormalised, cutoff-independent ghost, 
gluon and quark propagators. 

Together with the truncations for the other vertices, specified in appendix [Bj we end up with the following set of 
coupled, nonlinear integral equations for the ghost, gluon and quark propagators 

1 - Z 3 - g 2 N c [ dk K{p2 fj q2) G(k 2 )Z(q 2 ) , (11) 



G(p 2 ) ° a C J p 2 k 2 

1 _ „ , 2 N C /■ M(p 2 ,k 2 ,q 2 ) rin2 ^ f 2 . , 2 N C f ^^Q(p 2 ,k 2 ,q 2 )G(k 2 y- 2 -^G(q 2 )(- 2 -^ 



Zs + g 2 -^- / dfc ^ ; ' ' G(k 2 )G(q 2 )+g 2 -^ / dk 



Z(p 2 ) 3 J p 2 k 2 y ' vy ; y 3 J p 2 k 2 Z{k 2 f s Z{q 2 f s 

-*z* [ dk G ^ G ^ Mk 2 ) + Mq 2 ) 

g Z 2 2^ / M 2( k 2\„2 , M 2(„2\ 9 AJh2\AJn2\ l^J 



G{k) G(q) Mk 2 )+M<?) , 

k 2 + M 2 (k 2 )q 2 + M 2 {q 2 ) 2A l {k 2 )A l {q 2 ) 

R , 2l 7 L7 4 f ,, 47ra(g 2 )A M (g 2 ) 3g t (fc 2 ) 

Bi (p) = Z 2 m i + Z 2 - J dk k2A2{k2) + B2{k2y (13) 

A(2 , v , 7 4 f ^ ^a(q 2 )A u {q 2 ) A t (k 2 ) ( q 2 > P 2 + k 2 , (p 2 k 2 ) 2 \ 

Mp) = Z 2 + Z 2 - J dk - 2 ¥ ^ m — m ^-- + —^ + —^j , (14) 

where 5 = —9N C / (447V C — 8Nf) is the one-loop anomalous dimension of the ghost propagator and i S {u, d, s}. Note 
that with our choice l|lUfl of the quark-gluon-vertex it is the running coupling a(q 2 ) (cf. Eq.JSJl) multiplied with the 
quark self energy A u {q 2 ) that provides the interaction in the quark-SDEs. In order to allow the construction of a 
Bethe-Salpeter kernel in accordance with the axialvector Ward-Takahashi identity the interaction in the quark-SDEs 
has to be flavor independent and is therefore always chosen to contain the ^4-function of the lightest (up-)quark. The 
kernels ordered with respect to powers of q 2 = (p — k) 2 have the form: 



2\ 



1 / (p 2 -k 2 ) 2 \ 1 (p 2 + k 2 \ 1 



»r/ 2 ,2 2n 1 ( 1 2 fc2 lfc4 \ 1 lfc2 lf f 




= nnr\ - 



FIG. 2: Schwinger-Dyson equations for the quark and gluon propagator under the truncation scheme of the phenomenological 
model. All vertices are tree-level. Blobs represent fully dressed propagators. The seed term of the gluon equation is dressed 
with the Gaussian form Eq.El( see text). 



Q{p\k\q 2 ) = 



1 /lp 6 
7 \8k* 



9 - P V 



9p 2 A 9 k 2 
-— + 4 H 

4 k 2 ip 2 



W(p 2 ,k 2 ,q 2 ) = ^_+q 



1 

3p2 



1 k 6 
8 p z 



1 /V 2 2 

q z \ k A v 



1 1 

p2 y.2 

2fc 2 \ 2 k 2 k 4 
3p^J ~ 3 + V + V 



8p 2 fc 2 



(15) 



The momentum subtracted versions of Eqs. (|11|) - (|14[> are employed in our numerical investigations. (An outline of 
the subtraction procedure is given in section ITV A II ) 



B. The phenomenological model 

As a simpler (at least in some respects) alternative to the above described setup, we employ a model for the effective 
interaction between quarks and gluons, which has been developed in 0- This effective interaction can be related 
to the combination of dressed gluon propagator and dressed 7^-part of the quark-gluon vertex in the following way: 

T > a {p,k)Df v (p-k)T b v {k,p) — > iT a lfl g 2 Af u ( P -k)iT» 7 u (16) 

where both vertices are replaced by their tree-level counterparts {T°^ a — igT a 7il ). The central premise of phenomeno- 
logically constructing the effective interaction is that of integrated infrared strength - the interaction, when input 
into the quark self-energy integral, must be capable of generating sufficient amounts of dynamical chiral symmetry 
breaking [D^SB] in order to generate the experimentally observed hadronic mass spectra. One such construction is 
the quenched, Landau gauge form 

g 2 A% ab (p) = S ab V(p) Air 2 d £ exp , (17) 

where t^uip) is the transverse momentum projector. The parameters d and u> are fitted by comparing the resultant 
pseudoscalar meson mass and leptonic decay constant to experiment. The parameter uj, which sets the position of the 
maximum of the interaction, can be regarded as a length scale associated with the interaction, and d relates to the 
overall magnitude. It turns out that the meson masses (and the pseudoscalar and vector leptonic decay constants) 
are largely independent of uj in the range uj ~ 0.2 — 0.7GeV. As an example uj = 0.5GeV, d = 16GeV~ 2 gives realistic 
results for the pseudoscalar and vector channels 0. Clearly, by fitting to experiment, there is implicit information 
contained within such a phenomenological interaction about quark loops contained within the choice of parameters. 
However, the quark loops are at least from the viewpoint of truncated equations an additional degree of freedom. In 
the context of a full field theory (QCD here) such inclusions are self-consistently present as described in the previous 
subsection. In order to make such contributions also apparent in this phenomenological model we unquench the 
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effective interaction by adding to the phenomenological Gaussian form a quark loop term which generates an infinite 
series of increasing numbers of quark loops. The quark-gluon vertex itself remains undressed. The truncated set of 
SDEs for this phenomenological model is shown in Fig. [3 
The equation for the unquenchcd effective interaction is 



A2t(p) 



,0,oi 



(P) 



A ^(p)g 2 S cd t Xp (p)p 2 J:(p 2 )Af l/ (p) 



with 



] 5 cd t Xp (p)p 2 Z(p 2 )=g 2 S cd t Xp (p)p 2 Z A + [ dk J2 Tr[igT c lx S l {k)igT d lp S\k~p)} 

J A A „ 



(18) 



(19) 



i—u,d,s 



where the sum runs over the light quark flavors and the trace is over color and Dirac indices. The prefactors in the 
definition of £ arc such that £ is a dimensionless scalar function. The renormalisation constant Z A is independent of 
the external momentum scale and has been introduced to cancel the UV divergence of the quark loop. One can see 
that the integral of the right-hand side in Eq. 1)19(1 is transverse with respect to momentum p in exactly the same way 
as for the one-loop perturbative quark loop contribution to the gluon polarisation due to the presence of the tree-level 
vertices. Making the color and tensor structure of the effective interaction explicit with 



Af v (k) = 6 a % v (k)A(k 2 



(20) 



and projecting out the scalar part A(fc 2 ) from Eqs. I|18|) and i|19fl (see section IIV A II for more discussion on the 
projection) we arrive at 



with 



S(p 2 ) = Z A 



dk 



E 

i—u.d,- 



A(p 2 ) 



A<V) 



l + . 9 2 A°(p 2 )p 2 E(j5 2 



k 2 + M 2 (k 2 ) q 2 + M 2 (q 2 ) A{k 2 )A(q 2 ) 



•2 



p-k 2 n P'k 



(21) 



(22) 



where the momentum q = k — p. 

For the quark dressing functions A(p 2 ) and B(p 2 ) we obtain 



B(p 2 ) = J dkg 2 A(q 2 



3B{k 2 ) 



A{p2) = 1 + 3i dk92A{q2) k^A 2 (k 2 ) + B 2 (k 2 



k 2 A 2 {k 2 ) + B 2 {k 2 ) ' 

A(k 2 ) ( q 2 p 2 + k 2 (p 2 -k 2 ) 2 



2p 2 



2 q 2 p 2 



(23) 
(24) 



Note that the renormalisation constants of the quark equations above have been set to unity since the effective 
interaction is exponentially damped in the UV and the subsequent integrals are convergent. These last four equations 
are used for our numerical investigations of unquenching effects in the phenomenological model. Further technical 
details concerning our method to solve these equations in the complex momentum plane (necessary to solve the 
Bethe-Salpeter equation in Euclidean space) are given in appendix IO 



III. THE BETHE-SALPETER EQUATION 

The homogeneous Bethe-Salpeter equation [BSE] for flavor non-singlet quark-antiquark mesons can be written 

r %&P) = J dkK a0 . 5l (p,k;P)[S i (k + )r^(k;P)S j (k.)] ig (25) 

where K is the Bethe-Salpeter kernel, momenta k + = k + £P and fc_ = k + (£ — 1)P are such that the total 
momentum P = k + — fc_ and £ = [0, 1] is the momentum partitioning parameter reflecting the arbitrariness in the 
relative momenta of the quark-antiquark pair. The flavor content is expressed through the indices i, j and for flavor 
non-singlet mesons the kernel can be written independently of flavor content. Greek indices (a . . .) refer to color and 
Dirac structure. The BSE is a parametric eigenvalue equation with discrete solutions P 2 = —M 2 where M n is the 
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mass of the resonance. The lowest mass solution corresponds to the physical ground state. Since P 2 is negative, the 
momenta k± are necessarily complex in Euclidean space and so the quark propagator functions must be evaluated 
with complex argument. This leads to technical issues that are dealt with in appendix O 

One of the defining characteristics of the light pseudoscalar meson sector is the pattern of chiral symmetry breaking 
and the identification of pions as the Goldstone bosons of the symmetry breaking. In the Bethe-Salpeter formalism 
this is associated with observance of the axialvector Ward-Takahashi identity [AXWTI] [5j. The AXWTI serves to 
relate the quark self-energy and the Bethe-Salpeter kernel such that the pattern of chiral symmetry breaking in the 
pseudoscalar meson sector is guaranteed, independently of truncation or model details. In the flavor non-singlet 
channel and (for illustrative purposes) with equal mass quarks the AXWTI is 

- iP^l{p;P) = 5- 1 (p+) 7 5 + 75^ 1 (P-) - 2m r T 5 (p;P) (26) 

where we have introduced the (color-singlet) axialvector and pseudoscalar vertices (r^ and T 5 , respectively) and m r 
is the renormalised quark mass. Defining 

A(p; P) = -iP^l (p; P) + 2m r r 5 (p; P) (27) 

then the inhomogeneous Bethe-Salpeter integral equations for and T 5 can be combined to give (replacing the color 
and Dirac indices) 

A a/3 (p;P) = [iZ 2l5 f + Z i 2m rl5 ] ap + J dk K a0;Sy (p,k; P)[S(k+)A(k; P)S(k-)]^ s . (28) 

Using the AXWTI to rewrite A in terms of the quark propagator gives 

[5- 1 (p + )75 +755- 1 (p-)] ct/3 = [zZ 27 5f + Z i 2m rlb ] ap + J Sk K a ^{p,k; P) [ l5 S(k-) + S(k+) l5 ] yS (29) 

and by writing the Schwinger-Dyson equation for the quark as 

S- 1 (p) = -Z 2 %f + Z 4 m r - E / (p) (30) 
one can eliminate the inhomogeneous terms to get 



[-E / (p + ) 75 -75E / (p_)] a/3 = J dkK a ^{p,k;P)[ l5 S(k^ + S(k+)'y 5 ] lS , (31) 

thereby making the relationship between the quark self-energy and the Bethe-Salpeter kernel K manifest. 

For the truncation scheme given by the full-SDE set, in order to satisfy Eq. (|31|l we specify the truncated Bethe- 
Salpeter kernel to be: 

K ap .^(p, k; P) -> [r°-'(p +> k+)] a7 [T»> b (k-,p-j\ Sfl 5 a % v ( P - k)Z 2 ^ aiP -J^fr ~ k) , (32) 

where A u is the vector self energy of the up-quark regardless of the remaining flavor structure, which is left implicit. 
In the case of the phenomenological model, the kernel is nothing more than the standard dressed ladder truncation: 

K af) . Sy (p,k;P)-+ [r°;> + ,fc+)] Q7 [r°- 6 (fc_,p_)] 5/3 5«V(P-fc)A(p-fc). (33) 

Notice that the addition of the explicit quark loop in the expression for the interaction does not alter the structure 
of the kernel. Expanding out the various factors (retaining the flavor content), the Bethe-Salpeter equation for the 
two schemes are written as 

r«(p;P) = -|z a I Jfe 4?ra(P ~ *y ~ k) V(P - k) lfl S\k + )ri (fc; P)ffl(fc_) 7 , (34) 

for the full SDE scheme and 

P%P) = -| J dkg 2 A(p-k)t^(p~k)^S\k + )T^(k;P)S^k^ (35) 

for the phenomenological model. 
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Meson channels are specified by the quantum numbers J PC and flavor content. We assume isospin degeneracy 
{m c i = m u ) and consider only flavor non-singlet channels such that the flavor content of the meson can be generated 
in the BSE by supplying the appropriate quark propagators S l (henceforth we drop the flavor indices). The J PC 
quantum numbers are input into the BSE by specifying the Dirac and Lorentz decomposition of the bound-state 
vertices T to have the relevant transformation properties. We shall be considering pseudoscalar (0 and vector 
(1 ) mesons whose decompositions are: 

T PS (p; P) = 75 [T (p; P) - i?T x {p; P) - ipT 2 (p; P) - [fj] T 3 (p; P)] , (36) 
rl(p;P) = tf[tr (p;P) + yr 1 (p;P)-tr 2 (p;P)+i[fJ]T 3 (p;P)] 
+pl [T 2 {p;P) + 2ifT 3 { P] P)] 

+pl [r 4 (p; P) + iPT 5 (p; P) - ipTe(p; P) + [?,j>] T 7 (p; P)] , (37) 

where 7J and pj^ are transverse to P M (the on-shell vector meson is transverse to its total momentum P) . The I\ are 
scalar functions of p 2 and p-P (P 2 is on-shell for solutions of the BSE and therefore fixed). The charge conjugation 
properties of the decompositions above are such that when considering charge eigenstates, the Tj are either odd or 
even under the interchange p-P — > —p-P. By projecting the BSE onto the covariant basis it becomes a coupled set 
of nonlinear integral equations for the IV The technical procedure of solving these equations is outlined in section 

The solution of the pseudoscalar BSE can be used to obtain the pion leptonic decay constant, / w , which will prove 
useful in fitting parameters to observables. In order to calculate /„• one must firstly normalise the BS vertex function 
r. The normalisation condition is derived by demanding that the residue of the pole in the four-point quark-antiquark 
Green's function (from which the BSE is derived) be unity 0,0]. For £ = 1/2 it reads 

2P^3Tr d f dk !r(k, -P) dS{ \ + p P,2) T(k, P)S(k - P/2) + T(k, -P)S(k + P/2)T(k, P) dS{ \~ P/2) } 

J I C-f> O^p, ) p2 = _ M 2 

(38) 

where the trace is over Dirac matrices, both propagators describe it-quarks and the conjugate vertex function T is 
defined as 

T(p,-P) = CT T (-p,-P)C- 1 (39) 

with the charge conjugation matrix C = —"j-ilA- The leptonic decay constant characterising the pion coupling to the 
point axial field is subsequently given by |44j 

U = j^Tr d J dkT(k,-P)S(k + P/2) l5 pS(k- P/2), (40) 

where again the trace is over Dirac matrices. Eq. (H()[l can also be applied to the kaon by inserting the corresponding 
quark propagators (with the apropriate normalisation). 

IV. NUMERICAL PROCEDURE AND RESULTS 
A. Technical details 

1. Renormalisation and subtraction procedure 

To regularise the coupled system of ghost, gluon and quark SDEs we apply a MOM regularisation scheme p(il2lIl46T |. 
If we write the equations (|TTj) . i(T2|) . l(T3|) and (|T4l symbolically as 

' Z 3 + U ghost (p 2 ), (41) 



G(p 2 ) 
1 



Z[ji2] = Z 3 + n glue (p 2 ), (42) 

A(p 2 ) = Z 2 + Z 2 n A (p 2 ), (43) 
B(p 2 ) - Z 4 m R + Z 2 n B (p 2 ), (44) 
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where ma is the renormalised quark mass in the Lagrangian of the theory, this procedure yields 

~ 2 ) = Q^j + u ahost(p 2 )-Ug hos t{n 2 ), (45) 
1 1 



G(p>) G(n 

+ Hglue(p 2 )~Tlglue(p 2 ), (46) 



2\ 



Z(P 2 ) Z{pL 

A&T) = i-^nx^J + n^), (47) 

B(p 2 ) = B( P 2 ) + 1+ ^ 2) (n B (p 2 ) - H^ 2 )) , (48) 

where /i 2 is a suitable renormalisation point. The a priori unknown renormalisation constants Z%, Z3 and Z-i drop out 
and instead one specifies normalisation conditions for G(/i 2 ), Z(/i 2 ) and A(/i 2 ). For a given renormalised coupling g 
the Slavnov- Taylor identity 

l = Z 1 =Z g Z 3 Zl' 2 . (49) 

enforces that one cannot choose the normalisation conditions for G(/i 2 ) and Z(/i 2 ) separately, but only for the 
product Z(n 2 ) G 2 (/i 2 ), see Ref. |20j. Multiplicative renormalisability implies that the dressing functions for different 
choices of [i are trivially related to each other by factors of the renormalisation constants. Note that the coupling 
a(p 2 ) = a(fj, 2 )G 2 (p 2 ) Z (p 2 ) and the quark mass function M(p 2 ) = B(p 2 ) / A(p 2 ) are independent of the renormalisation 
point. 

A further complication is the appearance of spurious quadratic divergencies in the gluon-SDE, which are not 
eliminated by the MOM procedure. In principle, quadratic divergencies can be avoided by contracting the right hand 
side of the gluon-SDE with the Brown-Pennington projector R^uip) = 5^ — 4^1^ |47L Esj. However, using this 
technique one picks up spurious longitudinal contributions from the right hand side of the gluon-SDE, which affect 
the value of k in the infrared power laws of the ghost and gluon dressing function, Eqs. Q [TtI |24|. We therefore 
employ an alternative procedure, which is well known in the context of exact renormalisation group equations |3lLl49j |. 
After contracting the right hand side of the gluon-SDE with the transverse projector, leading to Eq. 1|12(1 . quadratic 
divergencies show up as terms proportional to A 2 /p 2 , where A is the UV-cutoff. In the infrared, such terms can be 
identified unambiguously by fitting the right hand side to the form a/p 2 + b/(p 2 ) 2K . The coefficient a then measures 
all contributions from quadratic divergencies and is subtracted from the right hand side at each step of the numerical 
iteration procedure. Together with the MOM-scheme described above this procedure ensures that all renormalised 
dressing functions are independent of the regularisation scale 2 . 

The situation is somewhat simpler for the phenomenological model, specified in section III Bl The projection of 
the quark loop contribution to the interaction Eq. Ijl9|l can be carried out using the Brown-Pennington projector to 
eliminate quadratic divergences without further complication. Also, the form of the effective interaction Eq. (|21|l with 
the Gaussian term A°(p 2 ) alters the picture. A°(p 2 ) (and hence A(p 2 )) vanishes exponentially in the UV, leaving the 
integrals for Eq. and Eq. (2H UV convergent. The remaining integral for S(p 2 ) does though require regularisation 
and renormalisation since the quark functions do reduce to their tree-level values in the UV, producing a logarithmic 
divergence. Also here we use a subtractive scheme and define the renormalised dressing function S(p 2 ) to be 

E(p 2 ) = S(p 2 ) - S( M 2 ) (50) 

such that E(p 2 = /i 2 ) = for some spacelike renormalisation scale /i. We then take as our interaction 

A(p 2 ) = _ . (51) 

V ' l+g 2 A°(p 2 )p 2 Y,(p 2 ) V ' 

In what follows, we drop the overbar on X and the subtraction is implicitly assumed though not written explicitly 
for notational convenience. The renormalisation scale /i is in principle a free parameter of the model. However, since 
we are studying light mesons with an associated hadronic scale of ~ lGeV we expect /j, to be in this region. It 
turns out that there is only a slight dependence of the overall results on /i, reflecting the approximate observance of 



2 In Ref. Il7i quadratic divergencies have been eliminated by a different, formally less rigorous procedure, which nevertheless leads to the 
same results as the present approach within numerical accuracy. 
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renormalisation scale invariance (c.f. appendix [Dj) . Throughout the main body of the paper we use the single value 
// = lGeV. 

There is one major distinction between the full SDE set and the phenomenological model and this concerns the 
role of the quark mass m^. In the full set of SDEs, the renormalisation of the quark mass and its subsequent running 
with momentum scales allows one to compare the current quark mass with the perturbative quark mass. In the case 
of the phenomenological model however, the quark mass is nothing more than a parameter of the model and cannot 
be directly compared to any other quantity. 



2. Solving the Bethe-Salpeter equation 

Projecting the BSE, Eq. |2Hlonto the covariant basis Ti (for pseudoscalar and vector mesons given via the decom- 
positions Eas. 1361 and 1371 respectively) it becomes a coupled set of nonlinear integral equations. To solve this set the 
Ti are approximated by a Chebyshev expansion in the angular variable z = p-P/ y/p 2 P 2 , 

N ch -1 

T l (p 2 ,z;P 2 ) = J2 ^ k T k {z)T k l {p 2 ■,P 2 ). (52) 

fc=0 

Since the T k form an orthonormal set one can further project the equations onto this basis. In the case of charge 
eigenstate mesons, one can extract an explicit factor p-P from the odd Ti and then use only the even T k . The meson 
mass solution to the system of coupled equations in the is then found in one of two equivalent ways. The first is 
to introduce a fictitious linear eigenvalue A(M) so that the original BSE reads 

\(M)T a p(p;P) = J dkK a0 . 5y (p,k;P)[S(k + )T(k;P)(k.)] 75 (53) 

with P 2 = —M 2 and vary the mass until X(M — M n ) — 1 (the physical ground state corresponds to the largest real 
eigenvalue). The second method is to recognise that after decomposition, the final integral over the radial momentum 
must be discretized giving rise to a matrix equation T^(p 2 \P 2 ) = K^f (pj , Pj, ; P 2 )T^, (p 2 ,;P 2 ) whose solution is the 
highest (physical ground state mass) P 2 = —M 2 for which Det(l — A") = 0. 

In principle, a truncation of the Chebyshev expansion destroys the Poincare covariance of the results, expressed 
through the invariance of the mass solution M 2 with the momentum sharing parameter £. However, the expansion 
quickly converges such that the numerical results are explicitly stable with varying £ 0. In this work we take £ = 1/2 
and N c h = 4(6) for the pseudoscalar (vector) meson channels. 



B. Numerical Results for the Propagator Functions 



Our results for the (un)quenched ghost and gluon dressing functions from the full set of SDEs are displayed in 
figure |3 In the infrared, the numerical SDE-results reproduce the analytical power laws, Eqs.JZJ. (This can be seen 
explicitly on a log-log-plot, displayed e.g. in Ref. H3)- I n the ultraviolet they reproduce the correct one-loop running 
of resummed perturbation theory, with anomalous dimensions 7 = 22n"-an s S an< ^ ^ = 44N 9 -8N f ^ or ^ e gl uon an d 
ghost dressing functions. Compared to the results of recent lattice calculations [33l I34I l36| we find good agreement 
for large and small momenta 3 . In the intermediate momentum region one clearly sees unquenching effects in the 
gluon dressing function. In this momentum region the system has enough energy to create quark- antiquark pairs from 
the vacuum. The screening effect from these pairs, first seen in Ref.|l8j|. drives the gluon dressing function closer to 
its perturbative form, i.e. the bump in the dressing functions becomes smaller. This effect is pronounced in both, 
the SDE and the lattice results. The overall difference in the size of the gluon bump between the two approaches 
indicates that the (neglected) gluon-two-loop diagrams in the gluon-SDE have an influence in this momentum region. 
Moreover, adding additional tensor structure to the quark-gluon vertex changes the magnitude of the unquenching 
effects slightly, as can be seen by comparison with Fig.l in Ref. [Hoj . In general, however, we expect that the inclusion 
of these contributions (in future work) will not change any of the conclusions drawn in this paper. 



Deviations for the value of the infrared exponent k may have systematical reasons, as discussed in l37l . These deviations are not 
important for the present investigation. 
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p [GeV] p [GeV] 



FIG. 3: Comparison of the quenched and unquenched ghost and gluon dressing functions with recent lattice data |33l l34l l36|. 
The sea-quark masses are m u / d ~ 16MeV, m a ~ 79MeV in the lattice simulations and m u / d ~ 3.9MeV, m a ~ 84MeV in the 
SDE-approach, c.f. table ITT1 



In comparing the quenched and unquenched interactions of the full SDE set and the phenomenological model, the 
situation is complicated by the presence of non-trivial quark masses. As pointed out previously, in the phenomenolog- 
ical model the quark masses are no more than parameters to be fitted to experimental results for meson observables 
(see next subsection for details). They cannot be compared directly to those of either the full SDE set or experiment. 
In addition, for both cases, the back reaction of the unquenching term induces changes to the fitted quark masses 
which further obscures the comparison. Therefore, to start the discussion we restrict for now to the chiral limit 
where a more unambiguous comparison can be made. Using the spacelike chiral quark propagator alone there are two 
quantities that can be compared in the quenched and unquenched truncations - the vacuum chiral quark condensate 
< qq >° and the quark mass function at zero momentum, M(0). The condensate is the order parameter for dynamical 
chiral symmetry breaking and at a renormalisation scale fi is defined as (see for example Ref. |5l| for a discussion on 
this topic) 

< qq >° (/i) = - lim Z 4 (fi,A)N c Tr d [ dkS {k) (54) 

where A is the regularisation scale, So is the chiral limit of the quark propagator and the trace is over Dirac matrices. 
In the full set of SDEs we determine the (un)qucnchcd condensate from the chiral quark propagator renormaliscd 
at a large scale /s 2 and subsequently evolve the value down to fj, = 2GeV employing one-loop scaling. For the 
phenomenological model, where Z± = 1, the condensate depends only indirectly on fi via the unquenching effects. 
In both approaches the unquenching effects on the quark-gluon interaction are calculated with physical sea quark 
configurations (see the next subsection for details of the fitting procedure and parameter values). 

Our results for < qq >° and M(0) are presented in tabled In the full SDE-setup with the quark-gluon vertex 
(I10|) we find a very slight enhancement of the condensate once quark-loops are taken into account. An even smaller 
effect has been found in Ref. 0]j where different approximations of the vertex have been investigated. These results 
together indicate that the condensate is almost constant for Nf < 3. Since the condensate is an order parameter for 
D^SB it is expected to change rapidly at the vicinity of the chiral phase transition. We conclude therefore that the 
critical number of flavors for the chiral transition of QCD at zero temperature is much larger than Nf = 3. This agrees 
well with the recently estimated value iVS = 10.0 ±0.4 in the framework of exact renormalisation group equations [&2 
and expectations from perturbation theory. For the quark mass at zero momentum we note a reduction of ~ 10% 
due to quark-loop effects in both the full SDE and the phenomenological frameworks. A similar reduction has been 
found in Ref. hj for a range of different truncations of the quark-gluon vertex. We therefore believe this effect to be 
model independent. Recent lattice simulations using staggered quarks confirm this result [53j. 

Next we compare the effective interaction of the phenomenological model, g 2 A(p 2 ) p 2 /4ir, to the interaction of the 
full-SDE setup, a(p 2 ) A(p 2 ), at finite quark masses. The parameters of the two schemes are fixed to reproduce the 
physical pion and kaon observables (see the next subsection for details of the fitting procedure and parameter values). 
From the plot in Fig. 01 one notes drastic differences. The interaction of the full-SDE setup stretches from the deep 
infrared to the ultraviolet, where the perturbative one- loop running is reproduced. In the infrared the interaction is 
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TABLE I: The chiral condensate and the chiral quark mass function at p 2 — in the full-SDE setup and in the phenomenological 
model. The unquenching effects are calculated with physical sea-quarks, see tab. [H] for the corresponding values of the quark 
masses and interaction parameters. 





-(< qq >°) 1/3 (MeV) 


M ch (p 2 = 0) (MeV) 


Model: 






N f = 


244 


366 


N f = 2 + 1 


243 


332 


full SDE: 






Nf =0 


266 


416 


N f = 2 + 1 


271 


388 




FIG. 4: Comparison of the quenched and unquenched interactions used in the two models. In the full-SDE framework the 
interaction is given by the product A(p 2 )a(p 2 ) whereas for the phenomenological model the equivalent function is g 2 A(p 2 )p 2 /4ty. 
The unquenching effects are calculated with physical sea-quarks, see tab. [2] for the corresponding values of the quark masses 
and interaction parameters. 

not vanishing due to the fixed point behaviour of the running coupling 4 , Eq. (JSJ, multiplied with a finite value of A(Q). 
The effective interaction of the phenomenological model, however, is localised at intermediate momenta. Despite the 
clear differences, both interactions produce similar qualitative effects in meson observables as will be seen in later 
sections. This is a validation of the assertion that the simplest qq meson states are dominated by the pattern of 
dynamical chiral symmetry breaking which arises from the integrated strength of the interaction and are not largely 
affected by the detailed structure of the interaction itself. 

The up/down and strange quark propagator functions for the full-SDE framework and for the phenomenological 
model are plotted in Fig. and Fig. HO respectively. Both sets of results for the quarks clearly exhibit a large degree of 
dynamical chiral symmetry breaking in the infrared region and we see the aforementioned effects of a reduction of the 
light quark mass function in the infrared once quark- loops are taken into account. For large momenta the numerical 
solutions of the full-SDE setup reproduce the logarithmic running known from resummed perturbation theory (see 
[Tc| for a detailed discussion of the ultraviolet properties of the quark mass function). Clearly, with the effective 
interaction of the phenomenological model being localised in the mid-momentum region, these asymptotic behaviours 
are not present. For intermediate momenta for all the quark flavors, we observe noticeable unquenching effects for 



4 Note that the infrared beh avio ur of the running coupling (defined from the Bjorken sum rule) has be en e xtracted very recently from 
experimental data at JLAB |. r )4l and found to be in excellent agreement with our SDE-result, taken from |17| . Since the exact theoretical 
relation between an invariant charge from experimental data and the coupling from the ghost-gluon vertex has not yet been clarified, 
such a comparison has to be treated with caution. Nevertheless, it suggests that the full SDE-interaction reflects essential properties of 
the strong interaction at large and small momenta. 
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FIG. 5: Comparison of the quenched and unquenched quark mass and wave-functions (left and right panels respectively) for 
the full-SDE setup. See text for details of the parameter configurations. 




FIG. 6: Comparison of the quenched and unquenched quark mass and wave-functions (left and right panels respectively) for 
the phenomenological model. See text for details of the parameter configurations. 



the full-SDE setup which are not present in the phenomenological model. In the infrared, the effect of unquenching 
is that Zf(p 2 ) is raised slightly and M(p 2 ) is lowered in all cases. There is, in addition, a striking difference in the 
behaviour of the quark wave function Zf(p 2 ) when comparing the two models (though the unquenching effects are 
similar to both). In the up-quark of the phenomenological model, the localised interaction results in bending Zf back 
upwards as one moves into the infrared region. This effect is suppressed for the strange quark. 

The last point of note concerns possible singularities in the quark functions in the complex plane. These can readily 
be identified with our technique for calculating the quark functions in the complex plane (see appendix Q , or using 
Schwinger- function techniques |55j . We find that the nearest singularities to the origin occur as pairs of complex 
conjugate poles in the timelike region. For the chiral quarks of the phenomenological model they are found at roughly 
(— 0.26 ± 0.17z)GeV 2 (quenched) and (— 0.23 ± 0.13z)GeV 2 (unquenched), whereas for the finite mass quarks they 
are further out in the timelike region. With the full-SDE setup the poles are further away in the timelike region as 
compared with the phenomenological model and there is little difference between the quenched and unquenched cases. 
The likely explanation for the full-SDE poles being further from the origin is again the observation that the associated 
interaction is not localised to a specific region but rather extends throughout the full momentum range. We conclude 
that the unquenching has a slight effect on the position of the singularities but will not affect the BSE calculations of 
later sections. Indeed it has been checked during the numerical calculations that no singularities appear within the 
boundaries of the integration region associated with the BSE for any of the results shown in the paper. 
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TABLE II: Parameter sets and results for m n , rn,K }k and m p for the quenched case (Nf = 0), the unquenched case with 
three degenerate 'sea'-quarks (Nf = 3) and the physical quark configuration case (Nf = 2 + 1) with two up/down quarks and 
one strange quark. In the phenomenological model uj = 0.5 GeV is held fixed and p = 1 GeV. In the full-SDE setup we give 
values for the quark masses read off at a large scale and evolved down to p = 2 GeV according to their one-loop behaviour. 





d(GeV" 2 ) 


m„(MeV) 


m s (MeV) 


m„(MeV) 


jW(MeV) 


m A -(MeV) 


/k (MeV) 


m p (MeV) 


model: 


















N f =0 


14.92 


5.26 


117.5 


139.6 


130.7 


493.7 


152.4 


737.6 


N f = 3 


16.94 


5.17 




139.6 


130.7 






734.7 


N f = 2 + 1 


16.66 


5.18 


116.4 


139.6 


130.7 


493.7 


153.0 


735.0 




AuoM(GeV) 


mt GeV (MeV) 


m 2 s GeV (MeV) 


m„(MeV) 


/W(MeV) 


m K (MeV) 


Jk (MeV) 


m p (MeV) 


full-SDE: 


















N f =0 


0.95 


4.17 


88.2 


139.7 


130.9 


494.5 


165.6 


708.0 


N f = 3 


1.16 


4.06 




139.7 


130.8 






690.0 


N f = 2 + 1 


1.09 


4.06 


86.0 


140.0 


131.1 


493.3 


169.5 


695.2 


PDG[56] 








139.6 


130.7 


493.7 


160.0 


770 



C. Numerical Results for Mesons: Physical quark masses, Nf = 2 + 1 

In this subsection we investigate two scenarios: quenched and unquenched with three physical quark masses, i.e. 
two light quarks representing the up/down quark and one heavier quark representing the strange quark. (Our results 
for three degenerate quark masses are presented in the next subsection.) 

To start, let us explain how we determine physical observables. In the case of the full-SDE setup we have to use one 
experimental quantity (here the pion decay constant f w = 130. 7 MeV) to convert the scale generated by dimensional 
transmutation to physical units. Furthermore we determine the quark masses m u / d and m s by the pseudoscalar 
observables — 139.6 MeV and mx — 493.7 MeV To roughly compare with the values of the particle data 

group we determine m u / d and m s from the quark mass function at large momenta and subsequently evolve them 

down to u = 2 GeV. To this end we extract the one-loop scale A N J^ = 0.95 GeV and A N JZ? = 1.09 GeV from 

MOM MOM 

the ultraviolet behaviour of the running coupling. Note that we have employed a MOM scheme, thus these scales 
have the expected magnitude 5 . In the case of the phenomenological model the scale is set by the paramater d in 
the effective interaction, where in addition we take the fixed value uj = 0.5 GeV and p = 1 GeV. In order not to 
obscure the present discussion, details of the procedure for fitting p are given in Appendix [D] The quark masses are 
determined by fitting to the experimental pion and kaon masses. 

The final values for the parameter sets with the resulting pseudoscalar meson masses, leptonic decay constants 
and vector meson masses are given in Table ITU When fitted to the experimental pion and kaon masses the resulting 
up/down and strange-quark masses are lowered when quark loop effects are taken into account. This has also been 
observed in corresponding lattice simulations [H^, EH . We furthermore see that the results for Jk and m p are quite 
insensitive to whether or not the system is unquenched. This leads to the conclusion that once the interaction has 
been fitted to the observable pseudoscalar observables, the vector meson mass is largely fixed. The likely explanation 
for this is that the ground state pseudoscalar and vector mesons are both states with the lowest orbital angular 
momentum (L = 0) in the sense of the naive (quantum mechanical) quark model - meaning that they are determined 
largely by the lowest spin contributions of the BS kernel given by the ladder approximation. The interaction then 
plays the same role in both channels, hence the similarity in results. Note that the p meson calculated within the 
framework of the truncated Bethe-Salpeter equation here refers to a pure quark-antiquark meson with no allowed 
decay channel - the physical p meson has a non-trivial decay width. Thus, one must allow for some change in the 
calculated mass if one were to take the decay process into account. However, our results for m p , especially in the 
full-SDE setup do seem low, although not unreasonable. 



5 For a discussion of the relation of a MOM-scheme to the AfS-scheme see section IV of Ref. l57l . Based on a three loop calculation the 

authors obtained the relation A-r-rrr — 0.346A . In our case this would result in the value A-t-ttt = 0.329 GeV for zero flavors. 

MS MOM MS 
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FIG. 7: Pseudoscalar (left panels) and vector (right panels) meson masses as functions of the quark mass parameter. We 
compare the quenched and degenerate unquenched cases for both the full-SDE (upper panels) and phenomenological model 
(lower panels). Recall that there is no direct comparison between the quark masses of the full-SDE and phenomenological 
model. 



D. Numerical Results for Mesons: Degenerate quark masses, Nf = 3 



Let us now consider the degenerate quark configuration case, with Nf — 3 identical quarks. To fit the parameters, 
we use a similar approach as that for the physical case, fitting to the pseudoscalar sector, though only the pion 
observables are relevant since there are no kaons. The parameters used and the resultant pion observables, along with 
the p meson mass are shown in Table [H] We see that the unquenched results for m p are only slightly changed by the 
different unquenching configuration. 

To further study the effects of unquenching we take the quenched and degenerate unquenched parameter configu- 
rations and plot the pseudoscalar and vector meson masses as a function of the quark mass parameter, keeping all 
other parameters fixed 6 . The results are shown in Fig. El We see that for the full-SDE setup, both the pseudoscalar 
and vector masses with larger quark masses are increased ~ 30MeV when we include the unquenching effects. For the 
phenomenological model we see the same trend though the effect is less pronounced. Notice that as m q — ► (the chiral 
limit) in both models, the pseudoscalar meson mass vanishes in accordance with the chiral symmetry considerations. 
This is a good check that the AXWTI has been faithfully maintained throughout the numerical procedure. 



The unquenched case with physical 'sea'-quark configurations is not considered because with the heavy meson masses generated by 
the 'valence'-quarks, the light quark propagators of the interaction need to be evaluated in a domain which includes the propagator 
singularities. This introduces technical problems that are beyond the scope of this work. 
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FIG. 8: Quark mass function at p — as a function of the quark mass parameter. We compare the quenched and degenerate 
unquenched cases. 



Given the differences in the models, one is led to investigate in more detail the relationship between the meson 
masses, the interaction and the quark mass functions and how this relationship is altered when the quark loop is 
included. We plot the quark mass function at p 2 — as a function of the quark mass parameter in Fig. El This serves 
to give a measure on the effect of unquenching for finite mass parameter quarks. We also plot in Fig. [5] the difference 
between the meson masses and twice the quark mass function M(p 2 = 0). In the static, quantum mechanical model 
of mesons this would be interpreted as a binding energy due to the interaction. We point out, however, that in the 
(relativistic) Bethe-Salpeter approach here this interpretation does not hold. We shall refer to this difference as the 
'binding' but the reader should bear in mind that this is a loose definition of the term. We see that in Fig. [5] the 
unquenching does indeed have a significant effect on the internal behaviour of the BSE. 

Expressed in the form of figs. 0D and one can now clearly see the difference in the two models. In the full-SDE 
setup at finite quark masses the quark mass function becomes insensitive to unquenching. This can be intuitively 
explained because the quark loops become suppressed by their mass. However, the overall meson masses show a 
~ 30MeV increase, reflected in the binding. Thus in the full-SDE setup, the unquenching of the interaction has an 
effect in the BSE but not in the quark mass function. 

In distinction, for the phcnomcnological model the situation is more complex. At low quark masses the quark mass 
functions and the binding do show the same behaviour as for the full-SDE setup - but of course for low quark masses 
all mesons are in the vicinity of the fixed physical point from which the parameters are fitted. As one increases the 
quark mass though, the curves for the quenched and unquenched quark mass functions and binding cross. Looking at 
the parameter sets, one sees that in the unquenched case the overall magnitude d is larger. As the quark masses are 
increased, the quark loop becomes suppressed and one is left with a stronger interaction, hence the increased quark 
mass function. (Recall that in the full-SDE setup there is no fitted parameter corresponding to the overall magnitude, 
hence the different behaviour of the quark mass function) . Since for small m q the unquenched mass function is lower 
than the quenched case eventually the curves cross. However, the meson mass does not increase so significantly under 
unquenching for larger m q , even though both the interaction and the quark mass function do increase. The resolution 
to this apparent dichotomy lies in the fact that for larger quark mass parameters, the relative contribution of the bare 
quark mass to the total meson mass is greater and the effect of increasing the interaction and quark mass function is 
damped (recall that the quark mass function referred to here is M(p 2 = 0)). 

This discussion of the detailed effects of unquenching serves to highlight that we are not dealing with a system where 
the total mass is given by the sum of the constituent masses and their binding energy but rather with a relativistic 
system whose constituents all depend non-trivially upon each other and crucially, all the components have varying 
effects at different momentum scales. The two models exhibit different behaviour when quark loops are included 
and this is because of the different forms of the interaction. The overall properties of the meson masses are largely 
determined by the integrated strength of the interaction. However, the closer inspection of unquenching for finite 
quark masses does indeed reveal that the form of the interaction plays a role. 

Finally, we make contact again to corresponding lattice calculations |58l l59| . Since it is most conclusive to compare 
'physical' quantities we plot the vector meson masses vs. the corresponding pseudoscalar meson masses. For large 
quark masses one may expect that unquenching effects in these quantities will be mainly determined by the interaction 
alone, since the decay of the vector meson to two pseudoscalar mesons is no longer kinematically allowed (though 
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FIG. 9: The difference between the pseudoscalar (left panel) and vector (right panel) meson masses and twice the quark mass 
function M(p 2 = 0) as functions of the quark mass parameter. We compare the quenched and degenerate unquenched cases. 



in principle unquenching effects involving vertex corrections beyond the truncation schemes discussed here may also 
play a role). The comparison to the lattice is therefore a genuine test of our two setups. From the plot in Fig. 1101 
we find good agreement between the curves from both of our setups and the available lattice data. It is apparent 
that the effects due to unquenching are smaller than both the differences due to the different truncation schemes and 
the systematic errors of the lattice results. Because we are plotting the more physical meson masses and not the 
(incomparable) quark mass parameters, the difference between the quenched and unquenched curves for both of our 
setups becomes the same - the curve for the unquenched vector meson mass as a function of the pseudoscalar meson 
mass is slightly raised for large values of the pseudoscalar meson mass. For small masses we observe a nonlinear 
relation of the two masses reflecting the different roles of the vector and the pseudoscalar meson sector in chiral limit. 



V. SUMMARY AND CONCLUSIONS 



We have investigated the effects of unquenching the gluon polarisation when applied to the Bethe-Salpeter frame- 
work for light pseudoscalar and vector meson masses. To obtain model-independent information on these effects we 
considered two distinct methods of calculating the corresponding gluon and quark propagators self-consistently. The 
first uses a sophisticated truncation of the Schwinger-Dyson equations for the ghost, gluon and quark propagator. 
The second uses a (simpler) phenomenologically motivated form for the quenched Yang-Mills sector of the theory. 
In both approaches, we explicitly took quark loop contributions to the gluon polarisation into account. Employing 
the resultant propagator functions in conjunction with a truncation for the Bethe-Salpeter equation we calculated 
light pseudoscalar and vector meson masses. The forms of the truncated Bethe-Salpeter kernel are such that the 
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N f =3, full SDE 




FIG. 10: Vector meson masses as a function of pseudoscalar meson masses. We compare the quenched and degenerate 
unquenched cases. The lattice results are taken from refs. [5sL l59| . 

axialvector Ward-Takahashi identity for color singlet, flavor non-singlet channels is obeyed, thereby ensuring that 
the generic features of the chiral limit of the pseudoscalar meson sector are manifest. Though the resulting effective 
quark-gluon interactions in our two approaches are rather different, both produce similar qualitative effects in meson 
observables. This finding enforces our general belief that the lowest meson states are dominated by the pattern of 
dynamical chiral symmetry breaking and are not largely affected by the detailed structure of the interaction. 

We observed sizeable unquenching effects due to the inclusion of quark loops to the gluon polarisation at various 
stages of our investigation. For the ghost and gluon propagator functions, calculated in the full Schwinger-Dyson 
truncation scheme, these effects occur predominantly in the intermediate momentum region. At momenta around 
lGeV the screening of quark-antiquark pairs moderates the bump in the gluon propagator considerably, whereas the 
infrared behaviour of both propagators is unaffected. For the quark propagator functions we found slight effects also 
in the infrared: the quark mass and wave functions for up/down and strange quark masses are reduced by about 
10% when unquenched. These results agree with previous findings an d have recently been confirmed in lattice 
simulations [33ll53j ]. We determined the location of the singularities of the quark propagator in the complex momentum 
plane and found a pair of complex conjugated poles, which are hardly affected by the inclusion of quark loop effects. 
Furthermore we found only very small changes in the chiral condensate as long as the number of flavors Nf < 3, 
indicating a large critical number of flavors for the chiral phase transition. All these effects are model independent in 
the sense that we observed them in the full SDE-setup as well as in our phenomenological model (where assessible). 

The main focus of the present work has been on the spectrum and decay constants of pseudoscalar and vector 
mesons. Here, the physical pion and kaon masses along with the pion leptonic decay constant have been used to 
determine the scale of the interaction as well as the up/down and strange quark masses. The resulting masses at 
the physical point are smaller whe n q uark loop effects are taken into account. Again this agrees with observations in 
corresponding lattice simulations 58, 59]. The kaon leptonic decay constant and vector meson mass were in rough 
agreement with available data, whereas the vector meson mass turned out slightly too low. The vector mass is subject 
to corrections due to a finite rho decay width |19| and further structure in the quark-gluon vertex |l5l ITfjl ] , which both 
were beyond the scope of the present work. Importantly however, our results for fx and m p were not significantly 
altered when quark loops were included. 

To further discuss the role of the quark loop contribution to the gluon polarisation, we compared quenched and 
unquenched mesons with degenerate sea-quark configurations. Keeping the scale of the interaction fixed, we studied 
the quark mass function as well as the pseudoscalar and vector meson masses as functions of the quark mass parameter. 
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As expected, both of our approaches preserve the Goldstone character of the pion in the chiral limit. For small quark 
masses the observables depend nonlinearly on the quark mass parameter. For large quark masses the two schemes 
display different behaviour under unquenching. The full Schwinger-Dyson truncation results exhibit an increase in the 
meson masses when quark loops are included, whereas for the phenomenological model this increase is much smaller. 

Furthermore we studied the vector meson masses as a function of the pseudoscalar meson mass and compared with 
recent lattice results from various groups. The resultant curves for the two schemes, though different from each other, 
are nonetheless in general agreement with the lattice data. For pion masses below 240 MeV, where no lattice data 
are available, our results show a nonlinear dependence of the vector meson mass on the pseudoscalar one. The effect 
of unquenching - when viewed as a function of the pseudoscalar meson mass - becomes the same for both schemes: 
the vector meson mass is slightly increased when quark loops are taken into account. This trend is also seen in the 
lattice simulations [58l l59j . where the effect is even more pronounced. However, these unquenching effects are small 
compared to the differences between both the truncation schemes we employed and the systematic errors of the lattice 
results. 
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APPENDIX A: CONVENTIONS FOR RENORMALISATION IN THE FULL SDE TRUNCATION 

SCHEME 

In the full Schwinger-Dyson equation truncation, renormalised dressing functions are denoted by an overbar. They 
are related to their unrenormalised counterpart via the renormalisation constant Zi(fi, A) where /i is the renormalisa- 
tion scale and A is the regularisation scale. The multiplicative renormalisation constants are defined as follows: 

S(p,A) = Z 2 (f,,A)S(p,fi) , (Al) 

D(p,A) = Z 3 (M,A)D(p,/i) , (A2) 

G(p,A) = Z 3 (M,A)G(p, M ) , (A3) 

T p ( Pl ,A) = z£(ji,A)r p (pi,n) , (A4) 

f„(p is A) = Z 1 - 1 0i,A)r p (p i ,^) , (A5) 
r p „A(Pi,A) = Z^ l {n,K)V pvX {jPi,ij) . (A6) 

The renormalised fermion mass function is obtained with the relation Z 2 M\ = Z^M V and the running quark mass 
M p is independent of the renormalisation scale /i. The coupling is renormalised according to 



5 2 (A) = ^U 2 (M) (A7) 
z 3 z 3 



and the Slavnov- Taylor identity, which ensures the universality of the renormalised coupling, is expressed as 

z 1 = z, = z^ (Ag) 

Z\ Z\ Zip 



APPENDIX B: A TRUNCATION OF THE GHOST AND GLUON SDES 

Here we shortly summarise the truncation for the ghost- and gluon-SDEs as introduced in refs. [HEIU- In 
order to close the Eqs. and l|12H we have to specify expressions for the ghost-gluon vertex, the three-gluon vertex 
and the four-gluon vertex. Our ansatz for the quark-gluon vertex has been discussed in the main body of this work, 
cf. Ea.lfHty. 
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FIG. 11: Ghost-gluon vertex DS equation. 



The crucial property of the ghost-gluon vertex 

T^(p,q)=T t ,{p, q )f abc (Bl) 

is that both, incoming and outgoing ghost momentum factorise in Landau gauge |30j| . Indeed, factorisation of the 
outgoing ghost momentum p is immediate and occurs in all linear covariant gauges. Factorisation of the incoming 
ghost momentum q can be seen easily from its Schwinger-Dyson equation, Fig. 1111 and the transversality of the 
gluon propagator D^ Vl which implies that 1^0^(1 — q) = q^D^(l — q). Thus, provided the gluon-ghost kernel in 
Fig. ^2 is well behaved in the infrared (i.e. not too singular), we obtain a bare vertex for q^ — > 0. Indeed at least 
one self-consistent infrared solution of the whole tower of SDE exists (2(| , in which the gluon-ghost kernel indeed is 
singular, but not singular enough to spoil the argument. If this is true, then a bare ghost-gluon vertex 

T ll (q,p)=iq ll (B2) 

provides a good approximation to the full vertex in the infrared and also in the ultraviolet region of momentum. 
Indeed, recent lattice calculations |6lj| and also a study in the SDE-approach |62| support this approximation. 

One knows from the self-consistent infrared solution (2|| , that the ghost-loop is the leading diagram of the gluon- 
SDE in the infrared. Although the three-gluon and four-gluon vertices are singular in this scenario, they are not 
singular enough to overcome the suppression induced by the gluon dressing functions. Thus effects from the four- 
gluon vertex on the gluon-propagator are negligible in the infrared. As they are also subleading in the ultraviolet and 
moreover are technically complicated to evaluate we neglect these effects from the start by choosing the four-gluon 
vertex to be zero. 

On the other hand, although the effects from the gluon-loop including the three-gluon vertex are also negligible 
in the infrared, they contribute to the one-loop running of the propagator in the ultraviolet and have considerable 
impact on the mid-momentum behaviour of the propagator. One therefore has to include a suitable ansatz for the 
three-gluon vertex. Such an ansatz has been constructed and investigated in refs. 17, 18, 60] and is given by 



T;l\(k,p, q) = r%Z°(k,p, q)H 3g (k lPl q) (B3) 
where T a v q) is the bare vertex. The dressing function H% g is given by 

with the momenta q and p running in the loop. Here the anomalous dimension of the ghost propagator, 8 = 
—9N C /(4:AN C — 8Nf), together with the vertex renormalization constant Z\ ensure the correct ultraviolet running of 
the vertex with momenta and renormalization scale. Furthermore this choice leads to cutoff independent ghost and 
gluon dressing functions in the continuum. 



APPENDIX C: SOLVING THE SDE IN THE COMPLEX PLANE 



In the (Euclidean space version of the) Bethe-Salpeter equation, one is looking for a solution P 2 — —M 2 which 
necessarily introduces complex momenta into the problem. The quark propagator must be evaluated at momenta 



7 Note that the color structure of the ghost and gluon propagators is just a trivial delta-function. Thus the color parts of all dressed 
Yang-Mills vertices get directly contracted with the antisymmetric structure constant f abc of their bare counterparts in the respective 
loops of the ghost and gluon SDEs. Possible color-symmetric contributions to the fully dressed vertices therefore do not contribute in 
these equations. We therefore use only the antisymmetric structure constants from the start. 
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within a parabola k 2 — M 2 /4 ± i\Jk 2 M 2 in the complex plane. In order to achieve this, we must solve the coupled 
quark and gluon propagator SDEs in the complex plane and we shall do this using a combination of contour integrals 
and UV asymptotic expansion. Expanding the integral measure and paying attention to the flavor index i, the coupled 
equations (|24I23I22[) become 



a* (a;) 



3(2tt)3 



£(*) 



pOO pi 

J dyy J dz\/ 1 - z 2 g 2 A(y)v l {x + y- 2^/xyz) 

V - 2^/xyz) [3] 

y 



3^L + 2z 2 



A pOO pi 

b l (x) = m l + ^-)3 J Q dyy J ^ dz yfl - z 2 g 2 A{y)s\ 
E 3(2tt) 3 



dyy I dzyl — z 2 v l (y)v t (x + y — 2^/xyz) 



(1 -4z 2 ) + 3 



y/v. 



(CI) 
(C2) 
(C3) 



where x — p 2 , y — k 2 and z = p-k/ y/p 2 k 2 . Alternatively, one can choose a symmetric momentum routing. Equations 
(IC1IC2IC3I) then become 



a l (x) 
b\x) 



3(2tt) 3 
4 



dyy dz\J\ - z 2 g 2 A(y-)v l (y+) 



xyz 



,/y x/A + yz 2 
y/x x/4 + y - yjxyz 



E 



3(2tt) 3 
4 



3(2tt) 3 



dyy I dzVl-z 2 g 2 A(y_y(y + )l3] 

y 



dyy I dz\J\ - z 2 v\y-)v\y + ) 



(1 - Az z ) + - 
x y ' 4 



(C4) 
(C5) 
(C6) 



where y± = x/A + y± ^fxyz. The reason for the two momentum routings is the following: At large values of the radial 
integral variable y, the asymmetric momentum routing has factors of either g 2 A(y) or v l {y) which vanish quickly for 
large y. The fidelity of the angular integration is thus less important even when x and y are both large (large x here 
refers to spacelike momenta). In the symmetric case at large x and y, the angular integration involves integrating over 
functions which vary considerably over the entire momentum range and so is numerically much more difficult. For 
this reason, although the two sets of equations are formally equivalent (ignoring for now the possible complications 
of the UV-cutoff), numerically we use the asymmetric routing at large x such that the angular integral can be more 
readily done. The use of the symmetric routing comes about in the context of the contour integrals below. 

Let us now describe how the equations are solved in the complex plane. At large $lx the functions take asymptotic 
forms 



X — > OG 
= I 



fc=l 

N 



a(x) 
b(x) 

S(x) '"=' J " er loga; 



Ebk_ 
x k ' 

k=l 



N 

E 

fc=i 



(C7) 
(C8) 
(C9) 



For a given set of coefficients one can extrapolate these functions at arbitrary x to the complex region, as long as %lx 
is large. Now consider the contour defined by the lines (shown in Fig. I12|) 

z-(t) — —xq + t — 2i^/xQt, t = [s, A], 

z\{r) = —x + A + 2i^/x \r, r — [— 1,1], 

z+{t) = —x + t + 2iy/x~ot, t=[X,e], 

z s (r) = — x + e — 2iy/x a er, r = [— 1,1]. 



(CIO) 



With these definitions the integral measure along the closed contour is: 



dz 



dt(l — jyio/t) + 2i\J xq\ 



dr 



dt(l + ly/xn/t) - 2i^ 



L 



dr. 



(Cll) 



23 




FIG. 12: Sketch of the closed contour in the complex z plane. 



For z\ within the contour the Cauchy integral formula holds: 

2n/(zi) = jdz^-, (C12) 

so given the function on the contour one can derive its value at some arbitrary interior point. Notice that numerically 
this is easy to implement as long as z\ is not close to the contour itself. 

Turning to the above coupled system of SDEs we supply a set of starting functions a 1 , b % and £ along the contour 
and along the spacelike axis, from which one can obtain the asymptotic expansion coefficients. From this we can derive 
the functions v % , s % and g 2 A using the expansion (for large real parts of the argument) or the contour integration 
(for small arguments) at those values needed to form the integrands of the coupled equations. With this we can then 
recalculate the functions on the left-hand side of the equations for external x again around the contour or on the 
spacelike axis. This gives us an iterative procedure to solve the system of equations. As explained above, we use the 
asymmetric momentum routing for large 9?x. However, for small 3?x we use the symmetric routing since for some 
point x on the contour we now need the integrand functions at points within the smaller parabola associated with 
x/4 and hence well inside the contour. This ensures that the Cauchy integral formula can be reliably implemented. 
It also means that if v l and s 1 are singular in the complex region, we can spot this by choosing the contour such that 
the singularity lies between the contour of external x and the smaller parabolic domain accessed by the integrands. 

We point out that within the iterative procedure there is no interpolation needed, except for the asymptotic region 
where the functions are well known and slowly varying. Certainly an interpolating scheme could not, even in principle, 
be used in regions where there are singularities present. In addition, the contour integrals are quick, simple to set 
up and are more accurate than an interpolation. The last point is the usefulness of the method in solving the BSE - 
given the converged solution to the SDEs along the spacelike axis and around the contour, one can then obtain the 
quark functions at any point needed within the parabolic BSE domain directly. 

The various parameters associated with the numerical implementation are as follows. The tip of the contour is given 
by the parameter e = O.lGeV 2 (avoiding the intcgrable square-root singularity to make the numerics easier), whereas 
the end is given by A = llGeV . 3tx > lOGeV 2 is 'large' and we use the asymptotic expansion. There a numerical 
limit to how large the overall contour parameter, xo, which defines the size of the complex region that the system 
accesses can be and this comes from the requirement that the contour integral be reliable and reasonably fast. In 
practice it is found that xo < 0.3GeV 2 works well. If quark propagators are needed outside this domain then one can 
use the SDE with symmetric routing to evaluate the functions, given the solution within. This gives an upper limit to 
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TABLE III: Parameter sets and calculated < qq >° and f° (chiral limit) using the phenomenological model, uj = 0.5GeV is 
held fixed. 





d (GeV- 2 ) 


H (GeV 2 ) 


-(<W>°) 1/3 (MeV) 


^ (MeV) 


quenched 


16.0 




251.2 


118.9 




13.0 


0.2 


251.4 


119.2 




15.0 


0.5 


251.8 


118.9 


unquenched 


16.2 


0.7 


251.4 


118.3 


(Nf = 3) 


18.0 


1.0 


251.2 


117.7 




24.0 


2.0 


251.0 


116.2 




45.5 


5.0 


250.9 


113.1 



TABLE IV: Parameter sets and resultant m^, f n , vtik fn and m p for the phenomenological model, ui = 0.5GeV is held fixed. 





M (GeV 2 ) 


d(GeV~ 2 ) 


m u (MeV) 


m s (MeV) 


m^(MeV) 


U (MeV) 


m K (MeV) 


/a- (MeV) 


m p (MeV) 


quenched 




14.92 


5.265 


117.55 


139.6 


130.7 


493.7 


152.4 


737.6 


unquenched 


1.0 


16.66 


5.185 


116.4 


139.6 


130.7 


493.7 


153.0 


735.0 


(N f = 2 + 1) 


0.7 


15.12 


5.22 


117.0 


139.6 


130.7 


493.7 


152.8 


736.2 




1.3 


18.67 


5.14 


5.14 


139.6 


130.7 






733.2 


unquenched 


1.0 


16.94 


5.17 


5.17 


139.6 


130.7 






734.7 


(Nf = 3) 


0.7 


15.15 


5.22 


5.22 


139.6 


130.7 






735.8 



the domain of applicability of the technique but this is due to the numerical implementation, not the technique itself. 



APPENDIX D: DISCUSSION OF THE RENORMALISATION POINT DEPENDENCE OF THE 

PHENOMENOLOGICAL MODEL 

In the phenomenological model, the renormalisation scale [i is a free parameter to be fitted. Since the model is not 
fully multiplicatively renormalisable, the results do in general depend on fi. We show here how /i is constrained. Let 
us begin by considering the behaviour of the chiral condensate defined by Eq. (|54H and additionally an approximation 
to the pion leptonic decay constant in the chiral limit 0] 

= i £ "yj^w - ¥ ayb 'y - h M • ( D1 ) 

It is known that this approximation underestimates /„. by around 10%. We calculate /° for a range of [i whilst using 
< qq >° to fit d with w = 0.5GeV fixed in order to constrain the range of /i. The results are shown in Table ITTT1 
We see that there is a mild decrease in f° as fi increases. For large values of \i we see that d must be increased 
dramatically in order to generate the same degree of dynamical symmetry breaking given by < qq > . Conversely, for 
small values of fj, (< 0.7GeV 2 ), d must be smaller than the quenched value. This is not unnatural if one considers that 
in the model the typical scale is ~ lGeV - renormalisation at a scale greater larger or smaller is not physical. From 
the lattice results displayed in Fig. [3] and the results of the full-SDE truncation shown also in Fig. 0] it is expected 
that the unquenching of the system reduces the gluon dressing by ~ 10% and this indicates that /i ~ l.OGeV 2 is a 
physical choice. 

We present parameter sets and meson mass results for the physical (Nf = 2 + 1) and degenerate (Nf = 3) 
unquenching scenarios with varying values of /i in Table IIVI From the insensitivity of the rho meson mass and the 
kaon leptonic decay constant, it is clear that fitting [i to either of these results is not a good procedure and the 
preceeding discussion of the interaction strength is more relevant to constrain /i. The near-invariance of the meson 
masses with respect to fi is the remnant of the renormalisation scale invariance that observables should adhere to and 
that we observe it here is a good sign that the phenomenological treatment retains the character of the full theory, 
insofar as these observables are concerned. 
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FIG. 13: Vector meson masses as a function of pseudoscalar meson masses in the phenomenological model for different values 
of the renormalisation scale fi. We compare the quenched and degenerate unquenched cases. 



Finally (and for the sake of completeness), using the quenched and degenerate [Nf = 3) unquenched quark config- 
uration scenarios (using the parameters of Table llVjl we plot the vector meson mass as a function of the pseudoscalar 
meson mass in Fig. ^5] for increasing quark mass parameter with different values of /x (cf. Fig. I10|) . Between all the 
cases it is seen that the unquenching is a small, though non-negligible effect. For fi = 0.7GeV 2 , where the parameter 
d (from Table ITV|) is almost unchanged from the quenched case the unquenched results are almost identical to their 
quenched counterpart. For larger values of fi the effect due to unquenching becomes larger. 
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